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We study effects of tunnel coupling on a pair of parallel disk-shaped Bose-Einstein condensates 
with the self-attractive intrinsic nonlincarity. Each condensate is trapped in a combination of 
in-plane and transverse harmonic-oscillator potentials. It is shown that, depending on the self- 
interaction strength and tunneling coupling, the ground state of the system exhibits a phase 
transition which links three configurations: a symmetric one with equal numbers of atoms in 
the coupled condensates, an asymmetric configuration with a population imbalance (a mani- 
festation of the macroscopic quantum self-trapping), and the collapsing state. A modification 
of the phase diagram of the system in the presence of vortices in the disk-shaped condensates 
is reported too. The study of dynamics around the stationary configurations reveals properties 
which strongly depend on the symmetry of the configuration. 



1. Introduction 

It has been predicted in many theoretical works [l[ that two dilute symmetric Bose- 
Einstein condensates (BECs), which are weakly coupled by the tunneling of atoms 
across the separating potential barrier, can give rise to the macroscopic quantum 
self-trapping. In particular, in the case of attractive inter-atomic interactions, the 
ground state of the system shows a transition from the symmetric configuration 
(the Josephson regime), characterized by equal numbers of atoms in the coupled 
condensates, to an asymmetric state (the self-trapping regime), characterized by an 
imbalance in the number of atoms [2, 0j- In the case of repulsive interactions, the 
self-trapping occurs not in the ground state , but rather in the first antisymmetric 
excited state. In the latter case, the self-trapping was demonstrated in experiments 
with the condensate of 87 Rb atoms 0] (for a review, see Ref. 0]). 

The self-trapping in the BEC loaded into the double-well potential is a mani- 
festation of the general effect of the spontaneously symmetry breaking (SSB) in 
nonlinear systems. As said above, asymmetric states trapped in symmetric poten- 
tials are generated by SSB bifurcations from obvious symmetric or antisymmetric 
states, in the media with the attractive or repulsive intrinsic nonlinearity, respec- 
tively (the SSB under the action of competing attractive (cubic) and repulsive 
(quintic) terms was studied too, featuring closed bifurcation loops @, @j). In terms 
of BEC and other macroscopic quantum systems, the SSB may also be realized as 
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a phase transition, which replaces the original symmetric ground state by a new 
asymmetric one, when the strength of the self-attractive nonlinearity exceeds a 
certain critical value. A transition of this type was actually predicted earlier in 
classical systems, viz., in a model of dual-core nonlinear optical fibers with the 
self- focusing Kerr nonlinearity In connection to the interpretation of the SSB 
as the phase transition, it may be identified as the transition of the first or second 
kind (alias subcritical or supercritical type of the SSB bifurcation), depending on 
the form of the nonlinearity, spatial dimension, and the presence or absence of an 
external periodic potential (an optical lattice) acting along the additional spatial 
dimension (if any) 0, 0] . 

Theoretical studies of the SSB in BECs were extended in various directions, espe- 
cially for matter-wave solitons. In particular, the symmetry breaking of the solitons 
was predicted in various two-dimensional (2D) settings [8(, including the sponta- 
neous breaking of the skew symmetry of solitons and localized vortices trapped 
in double-layer condensates with mutually orthogonal orientations of quasi-one- 
dimensional optical lattices induced in the two layers A different variety of 
the 2D geometry, which gives rise to its own mode of the SSB, is based on a sym- 
metric set of four potential wells . Self-trapping of asymmetric states was also 
predicted in condensates formed of dipolar atoms, which interact via long-range 
forces [12], and in the context of the nonlinear Schrodinger equation with a gen- 



eral local nonlinearity [13|]. The symmetry breaking is possible not only in linear 
potentials composed of two wells, but also in a similarly structured pseudopoten- 
tials, which are produced by a symmetric spatial modulation of the non-linearity 
coefficient, with two sharp maxima 14], 15]. 

Another generalization was developed for the SSB in two- [H] and three- 
component (spinor) [TtJ BEC mixtures, where the asymmetry of the density pro- 
files in the two wells comes along with a difference in distributions of the different 
species. As concerns multi-component systems, the analysis of the SSB was also 
extended to Bose- Fermi mixtures (l8l |. 

On the other hand, it is commonly known that the self-attraction in BEC may 
cause collapse of the condensate in the form of a "bosenova" (in which case three- 
body recombinations become important, in addition to the usual two-body colli- 
sions [H]). Therefore, the SSB in BEC trapped in double- well (dual-core) potential 
may compete with the collapse. Recently, we have determined 13, 2l[ the domain of 



parameters of such a symmetric dual-core system above which the collapse occurs. 
In particular, the competition between the SSB and the onset of the collapse in 
a pair of parallel cigar-shaped atomic condensates weakly coupled by tunneling of 
atoms was investigated in Ref. [13]. Further, in Ref. (HJ, the SSB and collapse were 
studied in a quasi-lD bosonic Josephson junction made by a double-well potential 
in the axial direction, and by a harmonic potential in the radial directions. 

In the present paper we consider a different setup, namely, a pair of parallel disk- 
shaped atomic condensates weakly-coupled by tunneling of atoms and confined by 
harmonic-oscillator potentials. This setup is ideal to analyze the interplay of the 
nonlinearity and tunnel coupling in the presence of vortices in both condensates [3] • 
In contrast to Ref. [3], which described this system by a pair of linearly-coupled 
2D Gross-Pitaevskii equations (GPEs) with the cubic nonlinearity, and actually 
presented the analysis of the SSB only below the collapse threshold, in this work 
we use the more accurate system of equations with the nonpolynomial nonlinearity, 



specific to the 2D geometry 22j, and the competition of the SSB with the onset 
of the collapse is one of main goals. After formulating the model in Section 2, we 
consider, in Section 3, the ground state of the system, which shows a phase transi- 
tion between three possible configurations: a symmetric one, with equal numbers 
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of atoms in the two coupled condensates, an asymmetric configuration with a pop- 
ulation imbalance (the macroscopic self-trapping, in the present setting), and the 
collapsing state. Then, we perform a similar analysis for localized states carrying 
vorticity in each core, which changes the phase diagram of the system. Finally, we 
study the dynamics of the two disk-shaped condensates around the stationary con- 
figurations. Starting from a symmetric configuration, we predict small-amplitude 
Josephson-like oscillations, with periodic transfer of the population imbalance from 
one core to the other. Starting from an asymmetric configuration, we find, instead, 
large-amplitude oscillations, which preserve the population imbalance. The paper 
is concluded by summary and discussion of open problems in Section 4. 



2. The model 

2.1. The dimensional reduction from 3D to 2D 

The starting point is the three-dimensional GPE for the mean-field wave function, 
i/j(r,i), which describes BEC in two parallel identical disk-shaped traps separated 
by a potential barrier: 

= I {-—V 2 + mu 2 z [(z - z ) 2 + {z + z ) 2 ]\ ij) 
at 2 l m J 

+ + (1) 

m 

where z is the coordinate transversal to the disks, 2zq is the separation between 
their centers along the z-direction, the two harmonic potentials with frequency 
uj z account for the transverse trapping of atoms in each disk, and W(x, y) is the 
potential acting in the disk plane (it is assumed to be identical for both disks). As 
usual, a s is the s-wave inter-atomic scattering length [23I ]. 

The first objective is to reduce Eq. ((TJ) to a system of linearly coupled equations 
for 2D wave functions pertaining to the separate disks, $1,2- To this end, we modify 
the approach developed for the system of two parallel quasi- ID "cigars" in Ref. 
(iol |. adopting a superposition of two single-disk ansatze: 



exp 



(z-z ) 2 \ $i(x,y,t) 



+ exp 



^ 2 m (x,y,t) 2 J y/ m (x,y,t) 
{z + z ) 2 \ $ 2 (x,y,t) 



2r] 2 (x,y,t) 2 J ^/r] 2 (x,y,t) 



(2) 



where 771 (x, y, t) and r]2(x, y, t) are the thicknesses of the two disks along the z axis, 
and the ID part of each wave function is normalized to unity. 
We proceed by substituting ansatz ((2J) into the Lagrangian corresponding to Eq. 



2 
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+muj 2 [(z - z ) 2 + (z + z ) 2 

Airh 2 
m 



+ 2W(x,yM 2 + ^^\^}:. (3) 
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The underlying assumption is that distance 2zq between the disks is essentially 
larger than the size a z = \/h/ (moj z ) of the transverse confinement in each of them, 
2zq S> a z . Due to this condition, the part of the Lagrangian, which accounts for 
the tunneling and is produced by the overlap of the two components of the wave 
function in ansatz ([2]), if substituted into Lagrangian ([3]), takes the following form: 

L T = -K J y)$Z(x, y) + ^(x, y)& 2 {x, y)] dxdy, 

where the effective coupling coefficient is defined as 



2 

K = Hu,Ae X p(-%). (4) 



In fact, the main contribution to the linear coupling (tunneling) comes from region 
z 2 < r) 2 2 around the midpoint between the disks. In that region, the transverse- 
confinement radius is determined by the ground-state wave function of the ID 
harmonic oscillator, which has characteristic length a z in the z direction. 

Finally, the effective dynamical equations for the two linearly coupled disks (n = 
1, 2) are written as 



•at** 



-Vl + W(x,y +5— 

2 n, 
-(- 



+t 1 -2 + n. 



$ n - K $ 3 -n , (5) 



where the scaled interaction strength is g = V27T7, with 

7 = 2Na s /a z , (6) 

the scaled linear coupling is 

k = K/(hto z ), (7) 

and the respective axial widths are determined by algebraic equations, 

ri = l+g\$ n \ 2 r] n , . (8) 

Notice that in Eqs. ([5]) and (jHJ) we have used scaled variables, viz., the length 
measured in units of a z = yWj (muj z ), time in units of a;" 1 , and energy in units of 

?UJ Z . 

Exact solutions to Eqs. © can be found by way of the Cardano formula, 



1 A* - 12 



1 Al -12 , ,„ ( Al - 12\' 1/2 

±2^P , (9) 
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where the upper and lower signs correspond, respectively, to g > and g < 0, and 

1 /3 

A n = (3/2) 1/3 (9 5 2 |$ n | 4 + v / 3^256 + 27(7 4 |$ n | 8 ) . (10) 
2.2. Properties of the 2D model 

Thus, we have reduced the initial 3D problem, based on Eq. ([T]), to the 2D problem 
for the set of wave functions of the BECs trapped in the two disks, which obey Eqs. 
© and dS|). To stabilize 2D solitons and vortices, we choose the in-plane potential 
as that of the 2D harmonic oscillator, 

W(x,y) = ^\ 2 (x 2 + y 2 ), (11) 

where A is the adimensional frequency of the planar confinement. The system 
conserves the total number of atoms in the two disks, i.e., N\{t) + AT 2 (t) = 2 (in 
the scaled units), where 

N n = J \$ n (x,y,t)\ 2 dxdy (n=l,2). 

Also conserved are the total energy and angular momentum. 

It is relevant to mention that the (effectively) two-dimensional self-attractive 
condensate, trapped in a periodic optical-lattice potential, features not only the 
collapse, when its norm exceeds the corresponding critical value, but also dereal- 
ization, when the norm falls below a certain threshold (the latter effects is also 
known in other dimensions) [3]. In the present setting, the derealization does not 
occur, as we consider the situation with the harmonic-oscillator potentials confining 
the condensates in all directions. 

Vortex-soliton solutions to Eqs. ([5]) are sought for as 

where r and 9 are the polar coordinates in the (x, y) plane, and S is the integer vor- 
ticity. In this way, Eqs. ([5]) can be reduced to coupled nonpolynomial Schrodinger 
equations in the radial direction, 



dr 2 r dr J ~ r r 2 



= 2 



TT^ + — 7T~ + -o + AV + 2^ 



-(- 

2 Us 



K 3 _ n . (12) 



Further, stationary states are then obtained by setting cj) n (r,t) = u n (r) e l/int , with 
real functions u n , which obey a system of stationary radial equations, 



1 r ( d 2 1 d \ S 2 - 2 



ar z r or J r z r? n 



+ o ( -2 + % 



while widths 771,2 are still determined by Eqs. (j5J), with l^i^) 2 replaced by 2 - 
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Figure 1. (Color online) Norms N\ and (solid and dashed lines) of the ground (zero-vorticity, S = 0) 
state trapped in the two parallel disks, in the course of the evolution in imaginary time. The nonlincarity 
coefficient is 7 = —0.4, while the coupling constant is re = 0.25 (upper panel) and re = 0.15 (lower panel). 
The adimensional notation corresponds to the length measured in units of a z = ^/h/ (mu) z ), time in units 
of uiz, and energy in units of hw z . Recall that 7 and re are defined by Eqs. and l[7}, respectively. 



As said above, the main objective of the work is to predict the SSB of the 
symmetric solitons, with u\ (r) = U2 (r), rj\ = 772 - In the system of two linearly 
coupled GPEs with the usual cubic nonlinearity, this problem was studied in Ref. 
0]. Here, we seek for stationary solutions by means of direct simulations of the 
time-dependent cylindrically-symmetric coupled 2D equations (|12p . using a finite- 
difference Crank-Nicholson algorithm in the imaginary time |25j. The initial con- 
ditions were taken as 

$ n (r, t = 0) = C n r s exp (-Ar 2 /2) , (13) 

where C n (n = 1,2) are normalization constants. Note that Eq. (I13p gives the exact 
quantum-mechanical wave function of the stationary vortex configuration in the 
absence of the nonlinearity (g = 0) and linear coupling (k = 0). In our numerical 
simulations we choose C\ 7^ C2, but with C\ taken very close to C2, to initiate the 
development of the symmetry breaking, if it possible. In particular, the norms of 
functions <&i(r, t = 0) and ^(j, t = 0) are taken as 1.01 and 0.99, respectively. 



3. Numerical results 



3.1. The ground state and vortices 

In Fig. [T] we plot the norms iVi and N2 (solid and dashed lines) of the two cou- 
pled condensates, with zero vorticity, S = 0, in the course of the evolution of in 
the imaginary time, by choosing 7 = —0.4 [recall 7 is defined in Eq. ([6]), 7 < 
corresponding to the attractive interatomic interactions], and slightly asymmetric 
initial conditions (i.e. slightly imbalanced populations). As shown in the figure, 
with K = 0.25 (the upper panel) the symmetry is restored during the time evolu- 
tion, while with k = 0.15 (the lower panel) the asymmetry is strongly enhanced 
towards a finite population imbalance. 
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Figure 2. (Color online) Stationary density profiles of the ground state (S = 0) in the two disks (solid 
and dashed lines), corresponding to the cases shown in Fig. {TJ. The upper panels: radial density profiles 
p(r) (a), and axial density profiles p(z), see Eq. I|14[l . (b), for 7 = —0.4 and k = 0.25. Lower panels: the 
radial profiles (c) and the axial profiles (d) for 7 = —0.4 and ft = 0.15. Units are the same as in Fig.[T] 



In the framework of the 2D description, the factorized ansatz ([2]) yields the time- 
dependent radial density profile, p n (r,t) = |$ n (r, t)| 2 , and its axial counterpart, 

p n (z,t) = 2^ [~ J^e-^Wpnfat)]* . (14) 
Jo Vn{r,t) 

In Fig. [2] the corresponding final (stationary) density profiles are displayed in the 
two disks, by solid and dashed lines. In the upper panels, the density profiles of the 
symmetric state are fully superimposed, while in the lower panels they are clearly 
distinguishable, for the asymmetric mode. 

Results of a systematic analysis, generated by varying parameters 7 and k, are 
summarized in Fig. [3j Here we show the phase diagram generated by the linearly- 
coupled system of 2D equations (fl2|) . with S = 0, in the parameter plane. As 
explained also in the the caption to the figure, in regions "symmetric" and "SSB" 
the system supports, respectively, stable symmetric and asymmetric stationary 
solutions. In region "collapse", the imaginary-time dynamics evolves towards a 
configuration with a zero- length axial width (in one or both disks). Notice that, on 
the right side of the vertical dashed line in Fig. [3] (i.e., at I7I > 1.07), the system 
always suffers the collapse, as in that region the nonlinearity strength exceeds the 
critical value leading to the onset of the collapse. 

The density profiles of the vortical states with S = 1 in both disks are displayed 
in Fig. [31 by means of the solid and dashed lines. The left panels of the figure 
clearly show the impact of the vorticity on the radial profiles p(r), which vanish at 
r — > 0. The figure also shows that, as expected, the transition to the asymmetric 
configuration follows reducing n. 

For the modes with 5 = 1, the phase diagram of the linearly-coupled system of 
equations (|12j) in the parameter plane of (7, k) is displayed in Fig. [5) Comparing 
Figs. [3] and [U we conclude that the "collapse region" is slightly reduced at the 
nonzero vorticity. In this case, the system always collapses at | —y | > 1.36. 



15 Molecular Physics 2dssb 



Taylor & Francis and I. T. Consultant 



0.8 



0.6 



SYMMETRIC 



*0.4 



0.2 




COLLAPSE 



1.2 



Figure 3. The phase diagram of the linearly-coupled system of 2D equations (1 1 L'[) with S = (the ground 
state). In regions "symmetric" and "SSB", the system supports, respectively, stable symmetric [«i(r) = 
U2(r)] and asymmetric [«i(r) ^ U2(r)] stationary solutions. The collapse takes place in the eponymous 
region. The system always suffers the collapse to the right of the vertical dashed line. Units are as in Fig. 
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Figure 4. (Color online) Stationary density profiles of vortices with S = 1 in the two disks (solid and 
dashed lines). Upper panels: (a) radial density profiles p(r) (a) and axial density profiles p(z) (b) for 
7 = —0.4 and re = 0.2. Lower panels: radial profiles (c) and axial profiles (d) for 7 = —0.05 and re = 0.15. 
The units are as in Fig. [T] 

The SSB can be characterized by the imbalance (asymmetry) parameter, 

N 1 (t)-N 2 {t) N 1 (t)-N 2 (t) 



N!(t) + N 2 (t) 



N 



(15) 



The competition between the symmetry breaking and collapse is further illustrated 
in Fig. [6] by plots of £ versus 7 for a relatively weak linear coupling, k = 0.1. In this 
figure, £ is the asymptotic value produced by the imaginary-time evolution in the 
framework of Eqs. (|12p with initial value £(0) = 0.01. The curves in Fig. [6] feature 
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Figure 5. The phase diagram of the linearly-coupled system of Eqs. (1 1 for the vortex modes with 5 = 1. 
In regions "symmetric" and "SSB", the system supports, respectively, stable symmetric [ui(r) = U2(r)} 
and asymmetric [ui(r) ^ U2(r)] stationary solutions. The collapse takes place in the eponymous region. 
The system always collapses on the right side of the vertical dashed line. Units are the same as in Fig. [T] 



a leap (represented by vertical segments) from the symmetric configuration with 
C = to the asymmetric one with £ ^ 0. Actually, the transition to asymmetric 
states in the present model always happens by a leap, i.e., the symmetry-breaking 
bifurcation is always subcritical, similar to the situation in the coupled equations 
with the self-attractive cubic nonlinearity 0]. The figure shows that, at fixed k, 
both the SSB and collapse happen at higher values of |7| in the case of S = 1, with 
respect to S = 0. 
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Figure 6. (Color online) The imbalance parameter, defined as per Eq. 1151 . as a function of interaction 
strength 7, for 5 = (filled triangles connected by the solid line), and 5 = 1 (open circles connected by 
the dashed line), for ft = 0.1. The curves terminate at the collapse points. Units are as in Fig. [T] 
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3.2. Real-time dynamics 

In the above subsections we have reported results of the imaginary-time simula- 
tions, which produce the stationary solutions. The next step is to test the stability 
of the modes by solving Eqs. (|12f) in real time. In Fig. [7] we display the real-time 
dynamics of the imbalance parameter, £(*)> °f the system with k = 0.1 for S = 0. 
The initial value is C(O) = 0.01. In the upper panel, we chose 7 = —0.1, which 
corresponds to a stationary symmetric configuration, while in the lower panel we 
set 7 = —0.4, which pertains to the asymmetric mode. 
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Figure 7. The real-time dynamics of the imbalance defined as per Eq. JTs}, for ft = 0.1 and S = 0. 
The upper and lower panels correspond, respectively, to values of the interaction strength 7 = —0.1 and 
7 = —0.4. The initial imbalance is £(0) = 0.01. Units are as in Fig. [T] 

Figure [7| shows that the dynamics are completely different in the two cases. With 
the initial imbalance £(0) = 0.01 in both cases, £(£) remains small in the course of 
the oscillations around the stationary symmetric configuration, changing its sign 
periodically. Actually, ((t) oscillates harmonically around the £ = 0. In the case of 
the stationary asymmetric configuration, the imbalance periodically assumes 
very large values, but it does not change the sign; actually, ((t) oscillates around a 
mean value, £ ~ 0.5. Note that the value of £ obtained asymptotically with these 
parameters (k = 0.1 and 7 = —0.4) in the imaginary-time simulations is £ = 0.89, 
see Fig. [6j 

In Fig. [8] we plot the evolution of £(t) for S = 1. It is important to stress 
that when the stability of the vortex solutions is tested against perturbations in 
real time, it is necessary to study the stability of the vortex against azimuthal 
perturbations, which may lead to splitting of the vortex that might seem stable in 
axially symmetric simulations [2ril |. For this reason, we employed full equations ([5]) 
to study the real-time dynamics of the vortices, considering both axially-symmetric 
initial conditions and those breaking the azimuthal symmetry. The general initial 
conditions used for the simulations of Eq. © are 



fcifo y,t = 0) = (x + iy)e- {x2+Sy2)/2 
& 2 (x, y,t = 0) = (x + iy)e- {x2+5y2)/2 



(16) 



We set 5 = 1 for the symmetric configurations, and 5 = 1.1 for ones with the broken 
azimuthal symmetry. Notice that $>i(x,y,t = 0) is normalized to N\(0) = 1.01 and 



15 Molecular Physics 2dssb 



Molecular Physics 11 



0.02 




n no ' 1 ' 1 ' 1 ' 1 ' 1 ' 1 ' 

' 50 100 150 200 250 300 
0.02 i 1 1 1 1 1 1 1 1 1 1 1 




n no ' 1 ' 1 ' 1 ' 1 ' 1 ' 1 ' 

' 50 100 150 200 250 30 
IF — 1 1 1 1 1 1 1 1 1 1 1 — =1 




Figure 8. (Color online) The real-time dynamics of the imbalance, defined by Eq. I] 1 5 p . for the vortex 
modes with S = 1, at K = 0.1. Four values of the interaction strength 7 are considered: (a) 7 = —0.1, (b) 
7 = —0.3, (c) 7 = —0.5, (d) 7 = —0.7. The initial imbalance is £(0) = 0.01, in all the cases. The solid lines 
are obtained by solving Eqs. <(5j for the axially-symmetric initial conditions l)16|l with 5 = 1, while the 
dashed lines are generated by the initial conditions with 8 = 1.1, which break the azimuthal symmetry. 
Units are as in Fig. [T] 



$ 2 (x,y,t = 0) to N 2 {0) = 0.99. In Fig. [SI the two upper panels [(a) with 7 = -0.1 
and (b) with 7 = 0.3] correspond to stationary symmetric configurations (the 
Josephson regime), with k = 0.1 (see Fig. [6|). The results displayed in these two 
panels of Fig. [8] show (the solid lines versus the dashed ones) that the additional 
azimuthal perturbation has no appreciable effects in the dynamics, apart from a 
slight dephasing. The third panel of Fig. [8] [(c), with 7 = —0.5] corresponds to 
a stationary asymmetric vortex (in the self-trapping regime). Here we find large- 
amplitude oscillations without a change in the sign of the population imbalance, 
C(t)- The solution with the unbroken azimuthal symmetry (the dashed line) has a 
period of oscillations very close to that observed in the solution with the azimuthal 
perturbation (the solid line). In any case, we conclude that the asymmetric vortex 
with k = 0.1, 7 = —0.5 and S = 1 is dynamically stable. Finally, in the bottom 
panel of Fig. [8]we plot £(t) for parameters S = 1, k = 0.1 and 7 = —0.7, which are 
at the border of the collapse region (see Fig. [6]). The panel shows that the solution 
is unstable (and eventually suffers the collapse). Here, the main difference between 
solid and dashed curves is the time after which displays the instability, which 
may be identified as the instant at which changes its sign for the first time. 
We observe that the instability with respect to the azimuthal perturbations can 
produce an additional border inside both the symmetric and asymmetric domains 
in Fig. [H We did not aim to produce this border in an exact form, as it is a 
computationally expensive objective. 
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4. Conclusions and open problems 



We have studied the dynamics of the self-attractive BEC in tunnel-coupled disk- 
shaped traps, by means of systematic simulations of the coupled nonpolynomial 
Schrodinger equations derived from the 3D Gross-Pitaevskii equation. In this way, 
we have investigated the phase diagram of the system as a function of the interac- 
tion strength (7) and tunneling coupling. We have found that borders of different 
domains in the phase diagram depend on vorticity S of the localized modes: both 
the SSB (spontaneous symmetry breaking) and collapse happen at larger values of 
7 in the case of S = 1 case with respect to the ground state (5 = 0). We have also 
studied the dynamics of the two disk-shaped condensates around the stationary 
configurations. Small-amplitude harmonic oscillations, showing a periodic transfer 
of atoms between the condensates, take place around the stable symmetric config- 
urations. Instead, large-amplitude oscillations without the change of the sign of the 
imbalance between the two condensates occur around the perturbed asymmetric 
configurations. 

There are many interesting open problems about Bose-Einstein condensates cou- 
pled by tunneling we want to face in the next future. In particular, we plan to inves- 
tigate quasi one-dimensional and quasi two-dimensional Bose-Einstein condensates 
in nonlinear lattices (i.e. with space-dependent interaction strength) [13] by us- 
ing the nonpolynomial Schrodinger equations. Moreover, we want to analyze the 
signatures of classical and quantum chaos [2^] in these double-well configurations. 
Finally, we aim to calculate analytically the coupling tunneling energy of bosons 
by means of the WKB semiclassical quantization [29J and comparing it with the 
numerical results of the Gross-Pitaevskii equation. 

LS thanks Luciano Reatto for 9 years of fruitful scientific collaboration at the 
Physics Department of the University of Milano. 
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